// EG4MagnetciField.cc -- implementation of the class
// created: Jan. 25, 2019, transcribed by Asia SUN, asia.rabbit@163.com
// Last modified: Jan. 25, 2019, by Aisa SUN

#include <fstream>
#include "EG4MagneticField.h"
#include "G4GenericMessenger.hh"
#include "G4SystemOfUnits.hh"
#include "G4UnitsTable.hh"

using std::cout;
using std::endl;

const double EG4MagneticField::kzMagIn = -555.5*mm;
const double EG4MagneticField::kzMagOut = 555.5*mm;

// trilinear interpolation
static double trilinear(const double*f, const double*dp);

EG4MagneticField::EG4MagneticField() : G4MagneticField(), fMessenger(nullptr),
fIsUniformMagneticField(false), fBy(0.*tesla), fIsFileLoaded(false){
  memset(fB, 0, sizeof(fB));
  if(!fIsUniformMagneticField) LoadMagneticFieldFile();

  // define commands for this clas
  DefineCommands();
} // end ctor

EG4MagneticField::~EG4MagneticField(){
  delete fMessenger;
} // end dtor

// position and time
void EG4MagneticField::GetFieldValue(const double p[4], double *bField) const{
  if(fIsUniformMagneticField){
    bField[0] = 0.; bField[2] = 0.;
    bField[1] = p[3] > kzMagIn && p[3] < kzMagOut ? fBy : 0.;
  } // end if uniform
  else GetMagneticIntensity(p, bField);
} // end member function GetFieldValue

void EG4MagneticField::DefineCommands(){
  // define /ETF/field command directory using generic messenger class
  fMessenger = new G4GenericMessenger(this, "/eg4/", "Field control");

  // field value command. Type: embedded class of G4
  auto& valueCmd = fMessenger->DeclareMethodWithUnit("field", "tesla",
    &EG4MagneticField::SetField, "set magnetic field strength of the dipole magnet");
  valueCmd.SetParameterName("field", true);
  valueCmd.SetDefaultValue("-1.");
} // end member function DefineCommands

// read the magnetic field measurement file, generated from measurements
// using Hall sensor, ASSIGN fB ARRAY
void EG4MagneticField::LoadMagneticFieldFile(const G4String &file){
  if(fIsFileLoaded) return;

  const char *filename = file.c_str();
  std::ifstream inFile(filename);
  if(!inFile.is_open()){
    cout << "ERROR: EG4MagneticField::LoadMagneticFieldFile():" << endl;
    cout << "File " << filename << " doesn't exist. Exiting...\n";
    exit(EXIT_FAILURE);
  } // end if
  G4String s_t;
  for(int i = 5; i--;) getline(inFile, s_t); // skip the commentary lines

  double current, x, y, z, Bx, By, Bz, Ihall, Tpt; // y(x)
  while(inFile >> current >> x >> y >> z >> Bx >> By >> Bz >> Ihall >> Tpt){
    if(fabs(Bx) > 1E50 || fabs(By) > 1E50 || fabs(By) > 1E50){
      cout << "EG4MagneticField::LoadMagneticFieldFile(): \
        Abnormal B value: (" << G4BestUnit(Bx, "Magnetic flux density")
        << ", " << G4BestUnit(By, "Magnetic flux density")
        << ", " << G4BestUnit(Bz, "Magnetic flux density") << ")"
        << endl;
    } // end if
    int i = G4int(-x/10.); // x: -1000, 0
    int j = G4int(y/10.); // y: 0, 90
    int k = G4int(z/10.); // z: -20, 1050
    if(i >= 0 && i <= 100
    && j >= 0 && j <= 9
    && k >= 0 && k <= 105){
      fB[0][i][j][k] = Bx;
      fB[1][i][j][k] = By;
      fB[2][i][j][k] = Bz;
    } // end if
    while(inFile.get() != '\n') continue; // skip the rest of the line
  } // end while
  if(inFile.eof()){
    cout << "EG4MagneticField::LoadMagneticFieldFile:"
      "End of the file reached.\nFile " << filename << " loaded successfully."
      << endl;
    fIsFileLoaded = true;
  } // end if
  else if(inFile.fail()){
    cout << "EG4MagneticField::LoadMagneticFieldFile\
      : Input terminated by data mismatch." << endl;
    fIsFileLoaded = false;
  } // end else if
  else{
    cout << "EG4MagneticField::LoadMagneticFieldFile: WTF?" << endl;
    while(1) getchar();
  } // end else
  inFile.close();
} // end of function LoadMagneticFieldFile

// assign B with the magnetic intensity vector at position (p[0], p[1], p[2])
// and time p[3]
bool EG4MagneticField::GetMagneticIntensity(const double p[4], double*B) const{
  // the position coordinates are deprived of their signs,
  // because of the symmetry of the magnetic field, the magnetic filed
  // intensity is irrelevant to the signs of the position coordinates.
  double pp[3] = {fabs(p[0]), fabs(p[1]), fabs(p[2])};

  // Expand the magnet field with measured B using extrapolation // ***********************
  // extrapolation for B at Y greater than 90mm and X greater than 390mm in the magnet
  // inside the cavity //
  if(pp[2] < 475.){
    if(pp[1] >= 90. && pp[1] <= 170.) pp[1] = 90. - 1.E2; // 135+35, including the V error
    if(pp[0] >= 390. && pp[0] <= 650.) pp[0] = 390. - 1.E2; // 620+30, including the H error
  }
  // between the iron wall (z=475) and the z=690 (coil) plane //
  if(pp[2] >= 475. && pp[2] <= 690.){
    if(pp[0] >= 390. && pp[0] < 700.){
      if(pp[2] - 690. < pp[0] - 390.) pp[2] = 690. + 1.E2; // z-ward extrapolation
      else pp[0] = 390. - 1.E2; // x-ward extrapolation
    }
    if(pp[1] >= 90. && pp[1] < 400.) pp[1] = 90. - 1.E2; // FEEL MY ECSTASY ~~~~~~~
  }
  // out of the inner area //
  if(pp[2] > 690.){
    if(pp[1] >= 90. && pp[1] < 500.) pp[1] = 90. - 1.E2; // AND YOU TOGETHER - ECSTASY ~~~
  }
  // Expand the magnet field with measured B using extrapolation // ***********************

  int ii[3] = {G4int(pp[0]), G4int(pp[1]), G4int(pp[2])};
  for(int &x : ii) x /= 10;
  double dp[3] = {pp[0]-ii[0]*10., pp[1]-ii[1]*10., pp[2]-ii[2]*10.}; // for interpolation
  double B_t[3][8]{}; // [xyz][8 vetices]
  bool isOut = false;
  static const double BMax = 1.2111*tesla;
  const double scale = std::fabs(fBy / BMax);

  if(ii[0] >= 0 && ii[0] < 100
  && ii[1] >= 0 && ii[1] < 9
  && ii[2] >= 0 && ii[2] < 105){
    for(int i = 0; i < 3; i++){ // loop over x, y and z
      // plane in the back
      B_t[i][0] = fB[i][ii[0]][ii[1]][ii[2]];
      B_t[i][1] = fB[i][ii[0]+1][ii[1]][ii[2]];
      B_t[i][2] = fB[i][ii[0]][ii[1]+1][ii[2]];
      B_t[i][3] = fB[i][ii[0]+1][ii[1]+1][ii[2]];
      // plane in the front
      B_t[i][4] = fB[i][ii[0]][ii[1]][ii[2]+1];
      B_t[i][5] = fB[i][ii[0]+1][ii[1]][ii[2]+1];
      B_t[i][6] = fB[i][ii[0]][ii[1]+1][ii[2]+1];
      B_t[i][7] = fB[i][ii[0]+1][ii[1]+1][ii[2]+1];
    } // end for over i
    // XXX: the maximum field is 1.2111 tesla
    B[0] = trilinear(B_t[0], dp)*gauss*scale;
    B[1] = trilinear(B_t[1], dp)*gauss*scale;
    B[2] = trilinear(B_t[2], dp)*gauss*scale;
  } // end if(...&&...&&...)
  else{
    B[0] = 0.; B[1] = 0.; B[2] = 0.;
    isOut = true;
  } // end else

  return isOut;
} // end of function GetMagneticIntensity


// trilinear interpolation
double trilinear(const double*f, const double*dp){
  double fx[4], fy[2], fz; // temporary variables
  // interpolation along x-axis
  for(int i = 0; i < 4; i++)
    fx[i] = (10.-dp[0])/10.*f[2*i] + dp[0]/10.*f[2*i+1];
  // interpolation along y-axis
  for(int i = 0; i < 2; i++)
    fy[i] = (10.-dp[1])/10.*fx[2*i] + dp[1]/10.*fx[2*i+1];
  // interpolation along z-axis
  fz = (10.-dp[2])/10.*fy[0] + dp[2]/10.*fy[1];
  
  return fz;
} // end of function trilinear interpolation
